arXiv: 1504.00514v2 [astro-ph.GA] 3 Aug 2015 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 4 August 2015 (MN IATeX style file v2.2) 


Simulating the formation of massive seed black holes in the 
early Universe. II: Impact of rate coefficient uncertainties 


Simon C. O. Glover 

Universitdt Heidelberg, Zentrum fiir Astronomie, Institut fiir Theoretische Astrophysik, 
Albert-Ueberle-Strafie 2, 69120 Heidelberg, Germany 


4 August 2015 


ABSTRACT 

We investigate how uncertainties in the chemical and cooling rate coefficients relevant 
for a metal-free gas influence our ability to determine the critical ultraviolet field 
strength required to suppress H2 cooling in high-redshift atomic cooling halos. The 
suppression of H2 cooling is a necessary prerequisite for the gas to undergo direct 
collapse and form an intermediate mass black hole. These black holes can then act as 
seeds for the growth of the supermassive black holes (SMBHs) observed at redshifts z ^ 
6. The viability of this model for SMBH formation depends on the critical ultraviolet 
field strength, Jcrit: if this is too large, then too few seeds will form to explain the 
observed number density of SMBHs. We show in this paper that there are five key 
chemical reactions whose rate coefficients are uncertain enough to significantly affect 
t^crit- The most important of these is the collisional ionization of hydrogen by collisions 
with other hydrogen atoms, as the rate for this process is very poorly constrained at 
the low energies relevant for direct collapse. The total uncertainty introduced into 
•/crit by this and the other four reactions could in the worst case approach a factor 
of five. We also show that the use of outdated or inappropriate values for the rates 
of some chemical reactions in previous studies of the direct collapse mechanism may 
have significantly affected the values of Jcrit determined by these studies. 

Key words: astrochemistry - hydrodynamics - methods: numerical - molecular 
processes - cosmology: theory - quasars: general 


1 INTRODUCTION 

Observations of quasars at redshifts z > Q demonstrate that 
supermassive black holes (SMBHs) with masses of order 
10® Mq had already managed to form by the time that the 
Universe was one billion years old. However, this is diffi¬ 
cult to understand if these SMBHs form from the seed black 
holes with masses M ~ 10 Mq that are produced as the end¬ 
points of the evolution of massive stars. The problem is one 
of timescales: even if we assume that accretion onto these 
seed black holes is very efficient, the time required for one 
of them to grow to a mass of 10® Mq can easily exceed the 
age of the Universe at 2 : > 6 ( Tanaka fc Haiman|2009 1. 

For this reason, an alternative model for the formation 
of the progenitors of high-redshift SMBHs has recently at¬ 
tracted considerable interest. In this model, known as the 
direct collapse model, gas in a protogalaxy with a virial 
temperature Tvir > 10^ K that is illuminated by an ex¬ 
tremely strong extragalactic radiation field collapses quasi- 
isothermally, maintaining its temperature in the range T ~ 
5000-10000 K through the combined effe cts of atomic h ydro- 
gen cooling and H“ bound-free cooling ( Omukai|2001 1. The 


high temperature keeps the Jeans mass in the collapsing gas 
high and inhibits fragmentation during the collapse (see e.g. 
Bromm fc Loeb||2003| [Begelman, Volonteri, fc ReS||2006| 


Shang, Bryan fc Haiman|2010 1. The precise outcome of the 

collapse remains uncertain, but the most likely outcomes are 
either the formation of a short-lived supermassive star or a 
quasi-star, a massive, optically-thick core whose centre has 
already collapsed to form a black hole ([Begelman, Volonteri, | 


fc Rees[2006 Begelm^|2010 Schleicher et al.||2013 (. fn ei¬ 

ther case, the end result is the formation of an intermediate 
mass black hole (IMBH) with a mass M ~ 10^ Mq or larger. 

The time required for this much larger seed black hole 
to grow by accretion into a 10® Mq SMBH is short enough 
that this mechanism can explain the existence of the ob¬ 
served high redshift SMBHs. In addition, the strength of 
the extragalactic radiation field required to suppress H 2 for¬ 
mation and the associated cooling in the collapsing gas is 
much larger than the mean value produced by early stellar 


populations (Haiman, Abel fc Rees 2000 Greif fc Bromm 


20061. This means that the formation of an IMBH by di¬ 
rect collapse will be a rare event, occurring only in a small 
number of protogalaxies located in highly clustered regions 
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where the local radiation field strength can be orders of mag¬ 


nitude above the mean value (Dijkstra et al. 2008 Ahn et al. 
2009 Agarwal et al.|[2012l. The viability of the direct col¬ 


lapse model depends upon just how rare this outcome is, 
and hence on how strong the radiation field needs to be in 
order to sufficiently suppress H 2 formation. If direct collapse 
occurs too rarely, then not enough IMBH progenitors will be 
produced to explain the observed number density of SMBHs 
at high redshift. 

Unfortunately, there is little agreement on the precise 
value of the radiation field strength required. This is usually 
quantified in terms of the specific intensity of the ultraviolet 
portion of the radiation field at the Lyman limit. Follow¬ 
ing Haiman, Abel & Rees (2000 1 , it is common to measure 

^ and to write it 


this in units of 10 ergs ^cm ^ 


Hz- 


as J 21 . The minimum value of J 21 for which H 2 cooling is 
sufficiently suppressed is then commonly denoted as Jcrit- 
Physically, the value of Jcrit depends strongly on the shape 


of the spectrum of the extragalactic radiation field (Shang, 
Bryan & Haiman 2010[ [Sugimura, Omukai fc Inou^ 201^ 


Latif et al.|2015r Agarwal fc Khochfar||2015| [Agarwal et al.| 


2015 ) and on the strength of any extragalactic X-ray back¬ 
ground or cosmic ray background (Inayoshi fc Omukai |2011 


Inayoshi & Tanaka 2015 Latif et al. 20151. However, our 


ability to determine Jcrit accurately from simulations is also 
affected by uncertainties arising from the way in which the 
chemistry of H 2 is modelled in these simulations. 

There are three main ways in which the chemical model 
adopted can affect Jcrit. First, its value in some circum¬ 
stances is highly sensitive to the way in which the effects 
of H 2 self-shielding are modelled (IWolcott-Green & Haiman 


2011 Sugimura, Omukai & Inouel2014|. Second, for reasons 


of computational efficiency, it is common to model the chem¬ 
istry of the gas using a highly simplified chemical network, 
but if the model is made too simple, then errors of up to a 


factor of a few can be introduced into Jcrit (Glover 2015; 


hereafter. Paper I). Finally, as we explore in this paper, Jcrit 
is also sensitive to the values adopted for the rate coefficients 
of some of the key chemical reactions occurring in the gas. 

Since most chemical rate coefficients are not known ab¬ 
solutely precisely, uncertainties in the input rate coefficients 
inevitably introduce an uncertainty into the value of Jcrit 
predicted by the simulations. In this study, we attempt to 
quantify the size of this uncertainty and to identify which 
chemical reactions make the largest contributions to it. In 
addition, as many previous numerical studies of the direct 
collapse model have adopted values for some of the chemi¬ 
cal rate coefficients that are now known to differ significantly 
from the actual values, we also explore how sensitive Jcrit is 
to these choices. 

The structure of our paper is as follows. In Section[^ we 
briefly present the numerical method we use in our study. In 
Section we quantify the sensitivity of Jcrit to the rate of 
each of the chemical reactions included in our model. Using 
this information, we identify a subset of reactions for which 
the sensitivity is large and critically discuss the uncertainties 
that exist in the rate coefficients for each of these reactions. 
In Section we carry out a similar analysis for the differ¬ 
ent cooling processes included in our model. In Section 
we combine our results and attempt to estimate the overall 
uncertainty in Jcrit due to uncertainties in the chemical and 
cooling rate coefficients. We also compare the results of our 


analysis with those of previous studies of the direct collapse 
model. Finally, we summarize the main results of our paper 
in Section m 


2 NUMERICAL METHOD 
2.1 The one-zone model 


We model the thermal and chemical evolution of gravita¬ 
tionally collapsing gas in halos illuminated by a strong ex¬ 
tragalactic radiation field using the same one-zone model 
as described in detail in Paper 1. We assume that the gas 
density evolves as 


^ ^ _p 
dt tff ’ 


( 1 ) 


where tg = (37r/32Gp)^^^ is the gravitational free-fall time. 
The thermal evolution of the gas is followed by solving the 
internal energy equation 

de _ p dp 


dt p 2 dt 


+ r-A, 


( 2 ) 


where e is the internal energy density, p = (7 — l)e is the 
gas pressure, F is the radiative and chemical heating rate 
per unit volume and A is the radiative and chemical cooling 
rate per unit volume. These are computed using the detailed 
atomic and molecular cooling function introduced in |Glover| 
|fc Savin] ( |2009| ) and updated in Paper 1. 

To model the chemical evolution of the gas, we use the 
conservative reduced network derived in Paper 1. This con¬ 
sists of 25 reactions involving 8 different chemical species: 
e“, H’^, H, H“, HJ, H 2 , He and HeH''". In Paper I we showed 
that the difference between the value of Jcrit derived using 
this network and that derived using a comprehensive chem¬ 
ical model of primordial gas is around 1 %, much smaller 
than the uncertainties considered in this paper. The full list 
of reactions included in our model is given in Table 

We model the effects of the external ultraviolet back¬ 
ground using two different spectral shapes: a 10® K black- 
body, cut oft at hv > 13.6 eV to account for absorption by 
neutral hydrogen, and a 10^ K black-body cut off at the same 
energy. For brevity, we refer to these two spectral shapes in 
the remainder of the paper as T5 and T4 spectra, respec¬ 
tively. The normalization of the radiation field is specified in 
terms of its specific intensity at the Lyman limit, measured 
in units of 10“^^ ergs“^cm“^Hz“^sr“'^, and written as J 21 . 

We assume that the gravitationally collapsing gas is op¬ 
tically thin in the continuum and model the effects of H 2 
self-shielding using the modified version of the |Draine fc 


Bertoldi (1996 1 shielding function given in Wolcott-Green, 


Haiman & Bryan (20111: 


/ah(AH2,T) = 


0.965 


-b - 


0.035 


(1 + x/hy-^ (l-ba;)0-5 
X exp [—8.5 X 10~"^(1 -|- a:)®'®] 


( 3 ) 


where x = Ahj/S x lO^'^ cm“^, 65 = 6 / 10 ® cm s“^, and 6 is 
the Doppler broadening parameter, which we assume to be 
dominated by the effects of thermal broadening. To compute 
Ahj , we assume that the dominant contribution comes from 
within a central dense core with radius equal to the Jeans 
length, and hence estimate Nh 2 as Aua = where Aj 

is the Jeans length. 
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As in Paper I, we caution that the simple approximation 
used here to compute A^H 2 may yield values that are system¬ 
atically incorrect by up to an order of magnitude. The abso¬ 
lute values derived here for Jcrit should therefore be treated 
with caution. However, this systematic uncertainty should 
not affect the conclusions we draw here regarding the sensi¬ 
tivity of Jcrit to uncertainties in the chemical rate coefficient 
adopted in the model. 


2.2 Initial conditions 


In paper I, we found, in common with previous authors, 
that the value of Jcrit depended only very weakly on our 
choices for the initial temperature, density and ionization 
state of the gas. The reason for this lack of sensitivity is 
that the question of whether or not enough H 2 forms in the 
gas to provide efficient cooling to temperatures T <C 10^ K is 
generally decided once the gas has already reached a number 
density of 100-1000 cm“^, by which time it retains only a 


weak memory of its initial conditions (see e.g. Shang, Bryan 
fc Haiman|2010 or Paper I). 

Therefore, in this study we consider runs performed us¬ 
ing only a single set of initial conditions. We set the initial 
temperature to To = 8000 K, the initial hydrogen nuclei 
number density to no = 0.3 cm“®, and assume that the ini¬ 
tial fractional ionization and H 2 fractional abundance are 
close to those found in the IGM prior to the onset of Pop. 
Ill star formation [xefi = 2 x 10“^, xh 2 ,o = 2 x 10“®). We 
assume that the gas is initially in charge balance, so that 
®H+,o = a^e.o, and place all of the hydrogen not accounted 
for by or H 2 in atomic form. We also assume that the 
helium starts entirely in neutral atomic form, and that the 
initial abundances of the other three chemical species in our 
network (H“, HJ, and HeH^) are zero. We note that these 
initial conditions are the same as those adopted in runs 2 
and 5 in paper I. 


3 IMPACT OF CHEMICAL RATE 

COEFFICIENT UNCERTAINTIES ON Jcrit 

3.1 How sensitive is Jcrit to the rate of each 
reaction? 

As the starting point for our study, we first explore how 
sensitively our derived value of Jcrit depends on the rate 
coefficient of each of the reactions in our chemical model. We 
do this for purely pragmatical reasons: if the value of Jcrit is 
insensitive to the value of a particular rate coefficient, there 
is little point in spending time and effort critically assessing 
how well we know that value of that rate coefficient. 

To assess how sensitive Jcrit is to variations in the rate 
coefficients of the reactions in our reduced network, we con¬ 
struct 52 different variants of our chemical model. In each of 
these variant models, the rate of one of the rate coefficients 
for our reduced set of 26 reactions is adjusted either upwards 
or downwards by a factor of 10®'® from its fiducial value. We 
then determine Jcrit for each of these models, for both a T4 
and a T5 spectrum. By comparing the outcome of the case 
in which a rate coefficient was adjusted upwards with the 
outcome of the case in which the same rate coefficient was 


adjusted downwards, we can quantify the effect of an order 
of magnitude change in that rate coefficient. 

The results of this analysis are shown in Table The 
sensitivity values we quote in the table are defined as the 
ratio of the largest and the smallest values of Jcrit that we 
obtain when modifying the rate coefficient of the specified 
reaction, i.e. 

c ..... Jcrit,max 

sensitivity = ——^- 

»>crit,min 

With this definition, in a case where Jcrit depends linearly 
on the rate coefficient in question we would expect to obtain 
a sensitivity of 10, while a value of 1 shows us that Jcrit is 
sensitive only to the inclusion of that reaction in the net¬ 
work, and not to the value of the rate coefficient adopted 
for the reaction. 

From Table we see that Jcrit is highly sensitive to the 
values of the rate coefficients for only a small subset of the 
reactions included in our minimal model. Specifically, if we 
restrict our attention to those reactions for which the sensi¬ 
tivity exceeds 1.25 (i.e. where an order of magnitude change 
in the rate coefficient results in a change in Jcrit of greater 
than 25%), then we find that in runs performed with a T4 
spectrum, the high sensitivity reactions are numbers 1-3, 5- 
8 , 10 and 12. In runs performed with a T5 spectrum, we 
find high sensitivities for a slightly different subset of reac¬ 
tions, numbers 1-3, 5, 6, 8, 10-12, and 22. This difference in 
behaviour is a consequence of the difference in the identity 
of the key process suppressing H 2 cooling in the runs with 
different incident spectra. In runs with a T4 spectrum, H“ 
photodetachment plays a much greater role than H 2 pho¬ 
todissociation in suppressing the H 2 abundance. Therefore, 
we find that in these runs, Jcrit is only weakly sensitive to 
the H 2 photodissociation rate (reaction 1), but has a much 
greater sensitivity to the H“ photodetachment rate (reac¬ 
tion 7). On the other hand, in runs with a T5 spectrum, H 2 
photodissociation dominates, and so the value of Jcut de¬ 
pends strongly on the rate of reaction 1, but has hardly any 
sensitivity at all to the rate of reaction 7, since in this case 
other reactions dominate the destruction of H“. 


3.2 Assessing uncertainties in the reaction rates 

In the previous section, we showed that Jcrit is highly sensi¬ 
tive to the values of the rates of only a small number of the 
reactions included in our reduced chemical model. In this 
subsection, we discuss how well we know the rates of these 
reactions. We restrict our attention to those reactions which 
have sensitivities greater than 1.25 for either the T4 or the 
T5 spectrum (or both). 


3.2.1 Photodissociation of H 2 (reaction 1) 

The molecular data on which the calculation of the H 2 pho¬ 
todissociation rate depends is known accurately (see e.g. 

), with an error of no more than a 
in principle it should be possible to 
compute the H 2 photodissociation rate produced by e.g. a 
T4 or T5 spectrum to a similar level of accuracy. However, 
in practice it is likely that the H 2 photodissociation rates 
used in studies of the direct collapse model are rather less 


Abgrall & Roueff 1989 


few percent. Therefore, 
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Table 1. Sensitivity of Jcrit to the rate coefficient adopted for each reaction in the reduced network 


No. 

Reaction 

Sensitivity (T4) 

Sensitivity (T5) 

1 

H2 -b 7 ^ H + H 

1.66 

9.52 

2 

H2-I-H-5-H-I-H-I-H 

16.0 

9.21 

3 

H- -b H -s. H2 -b e- 

6.31 

6.59 

4 

H+ -b H ^ H2 -b H+ 

1.10 

1.02 

5 

+ e“ ^ H + 7 

15.3 

18.1 

6 

H + e“ —>■ + 7 

7.01 

33.3 

7 

H“ + 7 ^ H + e“ 

5.54 

1.04 

8 

H-bH+ H+ -b7 

3.89 

1.40 

9 

H+ -b 7 -5- H+ -b H 

1.09 

1.00 

10 

H-bH-s>H+-be--bH 

1.34 

1.49 

11 

H--bH-s>H + H-be- 

1.12 

3.98 

12 

H -b e“ -s> H+ -b e“ -b e“ 

1.76 

1.67 

13 

H+ -b He -?• HeH+ -b H 

1.00 

1.00 

14 

H -b He H+ + e- + He 

1.10 

1.14 

15 

H2 -b H+ H+ + H 

1.00 

1.00 

16 

H2 -b He -5- H + H -b He 

1.00 

1.00 

17 

HeH+ -b H H+ -b He 

1.00 

1.00 

18 

H-bH-bH-i>H2-bH 

1.00 

1.00 

19 

H- -b He H -b He + e" 

1.00 

1.02 

20 

H+-bH-?-H-bH+-bH 

1.00 

1.00 

21 

He -b H+ -s. HeH+ -b 7 

1.00 

1.00 

22 

H- -b H+ -s. H-b H 

1.05 

1.75 

23 

H+ -b e- -5- H-b H 

1.01 

1.02 

24 

HeH+ -b e- -5- He -b H 

1.00 

1.00 

25 

H- -bH+ H+ -be- 

1.00 

1.00 

26 

H” -b e” —H -b e” -b e” 

1.00 

1.00 


The sensitivity quantifies the change in Jcrit that occurs when we change the value of the rate coefficient by a factor of ten. 


accurate than this, owing to the approximations made when 
deriving them. 


One potential source of error comes from the value 
adopted for the H 2 photodissociation rate in optically thin 
gas (i.e. the rate in the absence of self-shielding). To com¬ 
pute this, it is necessary to specify the level populations of 
our collection of H 2 molecules. Most studies of H 2 photodis¬ 
sociation in primordial gas assume that the gas is at low den¬ 
sity and hence that the effects of rotational and vibrational 
excitation can be ignored. In that case, the only energy lev¬ 
els of H 2 that are populated are the J = 0 para-hydrogen 
ground state and the J = 1 ortho-hydrogen ground state. 
Some calculations of the H 2 photodissociation rate further 


assume that only the J — 0 state is populated (e.g. Abel 


et al. 19971, but in practice, if only the lowest rotational 


levels are populated, then the H 2 photodissociation rate is 
insensitive to the ortho-to-para hydrogen ratio. 


In dense gas, however, it is no longer reasonable to as¬ 
sume that only the J = 0 and J — 1 levels are populated. 
Instead, the populations of the excited rotational and vi¬ 
brational levels will approach their LTE values. In LTE, the 
H 2 photodissociation rate increases as the gas temperature 
increases, owing to the increasing importance of photodisso¬ 
ciation from levels with u > 0. This is illustrated in Figurej^ 
where we plot the H 2 photodissociation rate in the LTE limit 
as a function of temperature, assuming either a T5 spectrum 
(solid line) or a T4 spectrum (dashed line). We see that in 
the case of a T5 spectrum, the photodissociatiou rate of 
hot, thermally populated H 2 is around 30% larger than the 
rate at low temperatures, which in this case is the same as 
in the low density limit. In the case of the T4 spectrum, 



10 100 1000 10000 
Temperature [K] 


Figure 1. Temperature dependence of the H2 photodissociation 
rate in the LTE limit. Values are shown for both a T4 spectrum 
(dashed line) and a T5 spectrum (solid line). In both cases, we 
assume that J21 = 1. 


the temperature dependence of the rate is much stronger: 
it increases by roughly a factor of seven as we increase the 
temperature from 10 K to 10^ K. 

In order to help us assess how much influence the tem¬ 
perature dependence of the optically thin H 2 photodissoci¬ 
ation rate is likely to have on Jcrit, we have examined the 
extreme case in which we adopt the LTE photodissociation 
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Bertoldi (1996 1 self-shielding function, 


rate throughout the collapse of the gas. We hnd that in this 
case, Jcrit is reduced to 12.3 if we use a T4 spectrum or 
1210 if we adopt a T5 spectrum. Therefore, even in this ex¬ 
treme case, Jcrit changes by no more than around 30%. In 
practice, the error introduced into Jcrit by the use of the 
low density limit for the optically thin H 2 photodissociation 
rate will probably be smaller than this. This is because it 
is the behaviour of the gas at densities n ~ 1000 cm“® that 
determines the value of Jcrit, and as this is less than the H 2 

critical density, ricrit ~ 10 "^ cm“^, the H 2 level populations /sh orig(A^H 2 , T) = 
will not yet have reached their LTE values. 

A second effect that can potentially become impor¬ 
tant when the photodissociation rate is extremely strong is 
the photodissociation of H 2 from highly excited vibrational 
levels. These levels can be populated by UV pumping of 
ground state H 2 (Shull 1978 Draine fc Bertoldi] 1996 1 or by 
so-called formation pumping, i.e. the tendency for newly- 
formed H 2 to be in a highly excited rotational and vibra- 


in Wolcott-Green, Haiman & Bryan (20111. This version of 


the self-shielding function is intended for use when one has 
rotationally hot H 2 , i.e. when the H 2 rotational level popula¬ 


tions are in LTE (although Wolcott-Green, Haiman & Bryan 


| 2011 | assume that the vibrational level populations continue 
to be negligible). On the other hand, the original Draine & 


0.965 


0.035 


[l + x/b^Y'^ (1-1-a;)®'® 

X exp [—8.5 X 10~^(1 -I- , (4) 


tional state (see e.g. Launay et al. 19911. The typical life¬ 


times of highly ex cited vibrational states of H 2 are short, o f 
the order of 10® s ( Wolniewicz, Simbotin fc Dalgarno|l998 1 , 
and so photodissociation of H 2 from these states is signifi¬ 
cant only when the photodissociation rate is extremely large. 


This issue was investigated by Shull (19781, who showed 


that UV excitation of these highly excited vibrational levels 
becomes significant only for Lyman-Werner photon fluxes 
F > 3 X 10~'^ photons cm“^ s“^ Hz~^. If we convert this 
constraint to a constraint on J 21 , we find that for a T5 spec¬ 
trum it becomes approximately J 21 > 5 x 10®, while for a 
T4 spectrum we have J 21 > 10®. 

These values of J21 are much larger than the values of 
Jcrit derived using our one-zone model. However, one-zone 
calculations tend to underestimate Jcrit compared to full 3D 
calculations, because they are unable to capture the com¬ 
plex temperature and density structure of the gas in real 
atomic cooling halos (see e.g. the comparison of one-zone 
and 3D results in Shang, Bryan & Haiman 20101. 3D sim¬ 


ulations performed using a T5 spectrum typically find that 


Jcrit ~ a fewx lO"' (Shang, Bryan & Haiman 


2010 


Latif et al. 


2015), and models that account for the positive feedback 


from soft X-rays also find similarly high values (Inayoshi & 
Tanaka|2015 I. At these values of J 21 , radiative de-excitation 
of highly excited H 2 molecules remains more effective than 
UV excitation, but the latter effect could potentially con¬ 
tribute at around the 10% level. One effect of this will be 
to slightly suppress the H 2 formation rate, as some fraction 
of the newly formed, vibrationally hot H 2 molecules will be 
photodissociated before they can radiatively decay to the 
vibrational ground state. Another effect will be to slightly 
enhance the total H 2 photodissociation rate ( Shull|1978 ). A 
detailed accounting for the effect of this on the final value 
of Jcrit is far outside the scope of the present work, but it 
seems likely that the effect will be of the order of 10 % or 
less. 

Lastly, but most importantly, the values of Jcrit that 
we derive in simulations performed using a T5 spectrum are 
highly sensitive to the treatment of H 2 self-shielding adopted 


in the simulation (Wolcott-Green, Haiman & Bryan 2011 


Sugimura, Omukai & Inoue 2014). For ease of comparison 


to other recent calculations, the simulations presented in 
this paper were performed using the modified version of the 


with X = VH 2/5 X 10^'* cm“^ and 65 = 6/10® cms“^ gives 
a better fit when the H 2 is rotationally cold, with signifi¬ 
cant populations only in the J — 0 and J = 1 rotational 
levels. Depending on which version of fsh one chooses, one 
recovers very different values for Jcrit in runs with a T5 spec¬ 


trum (Sugimura, Omukai & Inoue 2014 1 . For example, in our 


model, using /ah,orig in place of fsh increases J^rit from 1640 
to approximately 9000. 

This prompts the question of which version of the self¬ 
shielding function one should adopt for these calculations. 
Unfortunately, this is not an easy question to answer. As 
we have already discussed, at the density at which Jcrit is 
determined in these simple one-zone calculations, the H 2 
level populations have not yet reached LTE. Nevertheless, 
the density is high enough that collisional population of lev¬ 
els with J > 1 is starting to become important, and so we 
cannot simply assume that the H 2 molecules are rotationally 
cold. It is likely that reality lies somewhere between these 
two limiting cases, but how far in between is difficult to say. 
Clarifying this will require one to track the evolution of the 
H 2 rotational level populations during the collapse, a task 
which lies outside of the scope of this present study. 

Finally, in addition to this sensitivity to the choice of 
H 2 self-shielding function, the value of Jcrit is also sensitive 
to the method used to determine Nh 2 - Most one zone cal¬ 
culations adopt the simple approximation Vh 2 = uhj Lchar, 
where Lchar is some characteristic length scale, but different 
authors define this different length scale in different ways. 
There is general agreement that it should be related to the 
Jeans length, but no consensus on whether Lchar = Aj or 
Lchar = Aj/2. Moreover, it is not always entirely clear what 
definition is being used for Aj - different definitions in the 
literature can easily disagree by a factor of two, depending 
on whether one derives Aj from the dispersion relation for a 
ID plane wave or by considering the balance between ther¬ 
mal and gravitational energy in a uniform density spherical 
cloud. Of course, in reality none of these prescriptions is cor¬ 
rect - [Wolcott-Green, Haiman fc Bryan] ( |2011[ ) have shown 
that the Jeans length approximation can yield values of N 112 
and fsh that differ by an order of magnitude or more from 
the true values that one derives by post-processing 3D sim¬ 
ulations. To eliminate this source of error, it is necessary to 
abandon the use of simple approximations for Nh 2 , and to 
compute it self-consistently within a simulation using an al¬ 
gorithm such as TreeCol ( Clark, Glover fc Klessen|2012|). 
An attempt to do just this is reported on by |Hartwig et al.j 
(20151, but their results lie outside of the scope of the present 


Draine & Bertoldi (1996) H 2 self-shielding function given 


study. 
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6 S. C. O. Glover 


3.2.2 CoUisional dissociation of H 2 by H (reaction 2) 


in Figure 2 of Esposito & Capitelli (20091 suggests that the 


As discussed in e.g. Martin, Schwarz & Mandy (19961, there 


are two separate processes that contribute to the rate of 
this reaction. First, there is direct collisional dissociation, in¬ 
volving a transition from a bound state of the H 2 molecule 
directly into the continuum of classically unbound states. 
Second, there is dissociative tunneling, which involves the 
excitation of the H 2 molecule into a quasi-bound state - a 
state that has an internal energy larger than is required for 
dissociation, but which is separated from the continuum by 
a barrier in the effective potential. H 2 molecules in these 
highly excited quasi-bound states dissociate only if they can 
tunnel through the barrier before they have a chance to de- 
excite to a fully bound state. At temperatures T < 4500 K, 
this second process dominates the total collisional dissoci¬ 
ation rate ( Martin, Keogh fc Mand^|1998 l, and it remains 
important even at significantly higher temperatures. Some 
previous studies of the direct collapse model have accounted 
only for the direct collisional dissociation of H 2 , and not for 
the influence of dissociative tunneling (see e.g. |Shang, Bryan] 
fc Haiman|2010 1 . In paper I, we showed that the neglect of 


dissociative tunnelling can lead to one overestimating Jcrit 
by up to a factor of two; similar results have also been re¬ 


ported recently by Latif et al. (2014|. 


Even if one does account for both processes, uncertain¬ 
ties in the rate coefficients adopted will also introduce some 
uncertainty into Jcrit. In our simulations, we adopt rate co¬ 
efficients for these processes based on the master equation 


study of Martin, Schwarz & Mandy (19961. They present 


complicated fitting functions for both rate coefficients that 
account for their dependence on both the temperature and 
the density of the gas. The Martin et al. study itself relies 
on state-to-state collisional rate coefficients and collisional 
dissociation rates for H-H 2 collisions computed by|Mandy fc| 


Martin (19931 using the quasi-classical trajectories (QCT) 
method and the LSTH potential energy surface ( Liu||1973{ 
[Truhlar fc Horowitz[|1978{ [Siegbahn fc Liu|1978p . 

The use of the QCT method instead of a fully quan- 
tal method introduces some inaccuracy into the rate coeffi¬ 
cients. At low temperatures, very large differences are found 


between the two methods (see e.g. Sun fc Dalgarno 1994 
[Boothroyd et al.|1996[ ), but much better agreement is found 
at high temperatures, and so the QCT results should be 
reasonably reliable at the temperatures relevant for H 2 colli¬ 
sional dissociation. Similarly, differences between the LSTH 
potential energy surface (PES) and more recent parameteri- 
zations of the H 3 PES (e.g.|Varandas et al.|1987 Boothroyd 


|et al.|[T991[ |1996| [Mielke, Garrett fc Peterson| 26 b^ are im¬ 

portant at low temperatures (compare, e.g., the results pre¬ 
sented in |Sun fc Dalgarno||1994| |Lepp, Buch fc Dalgarno| 
|1995| and |Wrathmall fc Flower||20OT I, particularly for pure 
rotational transitions, but appear to become much less im¬ 


portant at high temperatures (Martin, Keogh fc Mandy 


1998). In the absence of a comprehensive treatment of H 2 


excitation by H along the lines of the Martin et al. study 
but using collisional rates computed with a different PES, it 
is difficult to state with certainty the total error introduced 
into the collisional dissociation rate by the use of an older 
version of the PES. However, the comparison between dis¬ 
sociation rate coefficients computed using the LSTH surface 


error must be relatively small at the temperatures of inter¬ 
est in this study. A conservative estimate for the error would 
be 25%. This introduces an error into Jcrit of roughly 35% 
in simulations with a T4 spectrum and 25% in simulations 
with a T5 spectrum. 

3.2.3 Associative detaehment of H~ with H (reaction 3) 

The rate of this reaction has recently been measured by 
Kreckel et al. (20101 for temperatures in the range 1 < 


T < 10 K. The systematic error in these measurements 
is approximately 25%. The results agree well with the low 


temperature measurements of Gerlich et al. (2012) and the 


calculations of Cfzek et al. (1998 1 , differing by amounts that 


are less than this systematic error. We have investigated the 
effect of this uncertainty in the rate of reaction 3 by scaling 
the rate coefficient both up and down by 25% and deter¬ 
mining what effect this has on our estimate of Jcrit. In runs 
with a T4 spectrum, we find that Jcrit = 15.1 when we de¬ 
crease the rate coefficient by 25% and Jcrit = 21.6 when we 
increase the rate coefficient by 25%. In runs with a T5 spec¬ 
trum, the corresponding values are Jcrit = 1350 and 1930. 
The remaining uncertainty in the rate of reaction 3 therefore 
introduces an uncertainty of around 40% into our estimates 
of Jcrit. 

3 . 2.4 Radiative recombination of RG (reaction 5) 

The rate of this reaction is known extremely accurately and 
high quality analytical fits are available that describe the 
temperature dependence of the rate coefficient. For example, 
the fit that we use in our study for the case B recombination 
rate coefficient - taken f rom|Hui fc Gnedin|(|1997[ ) and based 
on the data presented in Ferland et al. (19921 - has a quoted 
fitting error of only 0.7% over the temperature range 1 < 
T < 10® K. Therefore, although Jcrit is highly sensitive to 
the rate of this reaction, the high accuracy with which the 
rate is known means that in general, this reaction does not 
contribute significantly to the uncertainty in our estimates 
of Jcrit. 

That said, one must still take care to use the correct re¬ 
combination rate coefficient. Gravitationally collapsing gas 
within an atomic cooling halo generally has a low fractional 
ionization, and hence the mean free path for ionizing pho¬ 
tons is small. For this reason, the on-the-spot approximation 
generally applies, and when modelling the recombination of 
H^ one should therefore use the case B recombination rate 
coefficient. If the case A rate coefficient is used instead (as 
in e.g. the original | Abel et al.||1997| chemical model for pri¬ 
mordial gas, or the recent simulations by Latif et al.||2015 1, 
and the ionizing photons produced by direct recombination 
to the atomic hydrogen ground state are not explicitly ac¬ 
counted for, then the net effect is to increase the recom¬ 
bination rate by roughly 60% in the temperature range of 
interest. To assess the effect that this would have on Jcrit, we 
have re-run our simulations, using the case A recombination 


rate coefficient from Hui fc Gnedin (19971 in place of the 


and the BKMP2 surface of Boothroyd et al. (1996 1 shown 


case B value. With this change, we find that Jcrit = 10.1 for 
a T4 spectrum and Jcrit = 871 for a T5 spectrum, demon¬ 
strating that the use of the wrong rate coefficient introduces 
an 80-90% error into our estimate of Jcrit. 
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Simulating seed BH formation 7 


3.2.5 Radiative association of H and e (reaction 6) 


Calculations of this rate in the literature generally derive 
it by applying detailed balance to the inverse process, H“ 
photodissociation (reaction 7). Since the H“ photodissoci¬ 
ation cross-section is known accurately, one would expect 
that there should be little uncertainty in the resulting rate 
coefficient. However, when constructing an analytical fit to 
the calculated rate coefficient, it is necessary to choose some 
range of temperatures over which to fit. If the results of these 
fits are then extrapolated outside of their range of validity, 
then there is no guarantee that the resulting rate coefficients 
will remain accurate, and in this case, different fits may then 
yield very different results. This effect is seen quite clearly in 
Figure where we compare several different analytical fits 
given in the literature. We plot results using the fits given 
in |Hutchins| ( 1976|) (dotted line), Abel et al. (19971 (dashed 
line) and Stancil, Lepp &: Dalgarn^ 1998[ ) (dot-dashed line), 
in each case showing the ratio of the values given by these 
fits to those given by the fit in Galli fc Palla| (19981. 

We see from Figure that all four prescriptions agree 
well for gas temperatures in the range 10 < T < 2500 K. 
Within this temperature range, typical differences between 
different versions of the rate coefficient are of the order of 
10% or less. At temperatures T > 2500 K, however, the 


Hutchins (19761 fit starts to diverge strongly from the other 


rates. This behaviour is unsurprising, as the range of validity 
quoted by |Hutchins| for the fit is 100 <T < 2500 K. 

Similarly, at T > 6000 K, the other three versions of the 
rate coefficient start to differ significantly: the |AbeI et al.| 
119971 version increases more strongly with increasing tem- 
perature than the [Galli fc Palla| (|1998[) version, while the 


Stancil, Lepp & Dalgarno (1998 \ version of the rate coeffi¬ 
cient increases less strongly than the Galli fc Pallaj (19981 


version. This difference in behaviour is presumably due to 
some or all of the analytical fits reaching the end of their 
temperature range of validity at T ~ 6000 K, although it is 
difficult to be certain as none of the other papers state the 
temperature range over which their fits are supposed to be 
valid. 

We have investigated the effect that this difference in 
rates has on the value of Jcrit- For runs with a T4 spectrum, 
we find that Jcrit = 16.9,18.1,20.0 and 25.4 when we use 
the rate coefficient from 


Stancil, Lepp & Dalgarno 
[Abel et al.| ( |1997| ), |Galli fc Palla| ( |1998[ ) or |Hutchins 


(19981, 


(19761, 


repsectively. For runs with a T5 spectrum, the corresponding 
values are Jcrit = 1500,1630,1970 and 2840. 

We know from our own work in Paper I and from the 


work of previous authors (e.g. Shang, Bryan & Haiman 20101 


that the value of Jcrit is primarily determined by the be¬ 
haviour of the H 2 fraction in gas with a density n ~ 10®cm“® 
and a temperature T ~ 7500 K. In these conditions, the 


rate coefficient for reaction 6 given by Hutchins (19761 is 


not valid, and consequently it is no surprise that the value 
of Jcrit that we recover when we use this rate differs signifi¬ 
cantly from those recovered using any of the other analytical 


fits. We do not recommend use of the Hutchins (19761 rate 


in future studies of the direct collapse model. 

As far as the other three fits are concerned, there seems 
to be no particular reason to prefer one over the others at 
the present time. We can therefore conclude that the current 
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Figure 2. Comparison of different determinations of the rate 
coefficient for the radiative association of H and e“ (reaction 
6). We show results for three different analytical fits, taken from 
|Hutchins| ||1976[l (dotted line), [Abel et al.| | |T997[ | (dashed line), 
and [Stancil, Lepp fc Dalgarno| ( |1998[ l (dot-dashed line). In each 
case, we plot the ratio of the rate coefficient given by the analytic 
fits presented in these papers with that presented in the review 
of jGalli fc Palla| l |1998| |. 



uncertainty in this rate coefficient introduces an uncertainty 
of between 20-30% into the value of Jcrit. 

Finally, it should be noted that the uncertainty in this 
rate coefficient may also affect the later thermal evolution of 
the gas in halos where J21 > Jcrit, as in this case, H~ bound- 
free emission becomes one of the most important cooling 
processes once the density exceeds n ~ 10® cm“® (see e.g. 
Omukai|2001 Schleicher, Spaans fc Glover||2010 (. However, 


a detailed study of this issue is outside of the scope of the 
present paper. 


3.2.6 Photodissociation of H (reaction 7) 

Accurate cross-sections for this process are gi ven byjde Jong 
(19721 and [Wishart] ( |1979[ ). More recently, [Miyake et al. 

(20101 have computed revised cross-sections that account for 
the effect of the prominent resonances in the cross-section 
at photon energies hv > 11 eV. However, the effect of these 
resonances on the photodissociation rate is small: for a T4 
spectrum, the rate is enhanced by no more than a couple of 
percent, while for a T5 spectrum, the enhancement is around 
20%. We have investigated the effect of this enhancement on 
Jcrit, and find that for the T4 spectrum, it is reduced by a 
few percent, to Jcrit = 17.7, while for the T5 spectrum, it is 
not significantly affected. It should also be noted that this 
small chemical uncertainty is dwarfed by the much larger 
astrophysical uncertainty arising from the fact that both 
the T4 and T5 spectra are relatively crude approximations 


to the spectra of real high-redshift protogalaxies ( 

Sugimura, 

Omukai & Inoue|2014||Agarwal & Khochfar|2015| 

Latif et al. 

2015 Agarwal et al.|2015 

1 . 
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Figure 3. Comparison of different analytical fits for the rate 
coefficient for the radiative association of H and H"*" (reaction 8). 
The dotted line shows the ratio of the fit given in |Galli &: Pall^ 
l |1998[ | to that in|Latif et al.| ||2015[|, w hile the dashed line shows 
the ratio of the |Shapiro &: Kang|l|1987| | fit to the |Latif et al.| l [2015[ | 
fit. Comparison of the|Latif et al.| fit with the rates tabulated in 
[Ramaker &: Peek] ( |1976| l shows that it agrees to within a few 
percent at the temperatures of interest. 


3.2.7 Radiative association of H and iT*" (reaction 8) 


Rate coefficients for this reaction have been computed by 
both|Ramaker & Peek| ( [l976[ l and |Stancil, Babb fc Dalgarno| 
(19931, and agree to within 3%. Unfortunately, there is not 


such a good level of agreement between the different fits to 
the [Ramaker &: Peek| (|1976[) re sults that have been used in 
the literature. Shapiro & Kang (19871 introduced a fit of the 
form 


fcs.SK = 1.85 X 10 - 23 ^ 1-8 cm%“^ 


( 5 ) 


for T < 6700 K and 


fe.SK = 5.81 X 10 


(iloo) 


-0.6657 log{T/56200) 


cm^s ^(6) 


for T > 6700 K that was later adopted by many other au¬ 
thors; for instance, it is this prescription that is currently 
used in the enzo primordial chemistry model, as used in e.g. 


Shang, Bryan & Haiman (20101. However, in their review of 


primordial chemistry, Galli & Palla (19981 use a different fit, 
ostensibly to the same data 

^2 


fcg.GP = dex [-19.38 - 1.523logT+1.118(logT)^ 


-0.1269(logT)®l cmS“\ 


( 7 ) 


Finally, in a recent paper, Latif et al. ( 2015| ) give the fit 

fcg.L = dex [-18.20 - 3.1941ogT+ 1.786(logT)^ 

-0.2072(logT)^] cmS"\ (8) 

citing Coppola et al. |( |MTT| as the source. 

In Figure we compare these three different versions 
of the rate coefficient by plotting the ratios feg,sK/feg,L and 
fcg,Gp/fc8,L as a function of temperature. We see that at the 
temperatures of interest, the scatter in the values of the 


three fits is around 30%, ten times larger than the difference 


between the Ramaker & Peek (19761 and Stancil, Babb & 


Dalgarno (1993) calculations. This scatter has an apprecia¬ 


ble effect on the values of Jcrit that we recover for runs 
performed using a T4 spectrum. We find that for runs using 


the Shapiro & Kang (1987), Galli & Palla (1998) and Latif 


et al. (2015) fits, Jcut = 19.8,16.8, and 18.0, respectively. 


a total difference of around 20%. A similar analysis carried 
out for runs using a T5 spectrum finds a much smaller range 
of values: for the same three cases, Jcrit = 1670,1600, and 
1630, giving us a total uncertainty in this case of only 4%. 

We have compared the values predicted by these three 
different versions of the rate coefficient with the values 


tabulated by Ramaker & Peek (1976), and find that the 


Latif et al. (20151 fit gives the best overall match to the 


tabulated data, agreeing to within around 4-5% over the 
entire temperature range considered by [Ramaker fc Pe^ 
(10 < T < 32000 K). We therefore recommend that this fit 
to the rate coefficient should be used in future numerical 
studies of the direct collapse model. 


3.2.8 Collisional ionization of H by H (reaction 10) 

Although there have been a number of calculations and 
measurements of the cross-section for this process at high 


collision energies (see e.g. McGlure 


1968 


Hill, Geddes fc 


Gilbody[[1979[ [Shingal, Bransden fc Flower[[1989| [Riley H 


Burke|1999 ), only a few studies have looked at its behaviour 


at the low energies relevant for the direct collapse model. 
In our fiducial model, we use the rate coefficient for this 


reaction given in Lenzuni, Ghernoff & Salpeter (1991): 


fcio.LGS = 1-2 X 10 ^^T^'^exp I 


157800\ 
T ) 


3 -1 
cm s 


( 9 ) 


This is based on the experimental cross-sections of [Gealy fc[ 
van Zyl (1987). However, their measurements only exist for 


collision energies E > 100 eV. Therefore, in order to derive 
a rate coefficient that is valid at temperatures T ~ 10^ K, 
it is necessary to extrapolate the measured cross-sections 
down to the ionization threshold energy of 13.6 eV. This 
extrapolation can potentially introduce a substantial error 
into the resulting rate coefficient. Indeed, a good example of 


this is provided by the fact that Hollenbach & McKee (19891, 


using the same cross-section data, derive a rate coefficient 

149000 \ 


fcio,HMS9 = 1.7 X 10-'^T° ®exp (- 


cmS“\ (10) 


T ) 

Although based on the same experimental data, these two 
expressions differ by a factor of 5 at T = 10^ K and by an 
order of magnitude at T = 7000 K. 

The H-H collisional ionization rate coefficient can 
also be derived using the semi-empirical method described 
in Drawin (1968 1969). Following the interpretation of 


Drawin s results given in Soon (19921, we can write the rate 


coefficient as 

fcio,D69 = 1.57 X cm^ s"\ 

where 




1 -I- (2me/mHa;)3 


exp(-a;). 


( 11 ) 


( 12 ) 


me is the electron mass, mn is the mass of a hydrogen atom, 
and a; = 157800/T. 
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Simulating seed BH formation 9 


Yet another expression for the H-H collisional ionization 


rate coefficient can be found in Hollenbach & McKee 119791. 


They assume that it is a co nstant factor of 1.7 x 10 
than the H-e“ rate (citing 


Drawin 


1969 


smaller 
in support of this 


fact), and hence write the rate coefficient as 

fcio,HM 79 = 9.8 X exp cm^ s“\ (13) 

Finally, analytical expressions for the H-H ionization 
cross-section and corresponding thermal rate coefficient are 
given in 1 


Soon (1992 l. These result from an extension of the 


work of Kune & Soon (19911, and were calculated using 


the classical impulse approximation (Gryzinski 19651. Un¬ 
fortunately, there are a couple of significant problems with 
the expressions presented in |Soon| s paper, as already noted 
by Barklem (20071. Most importantly. Soon uses the wrong 


value for the ionization threshold energy. In the centre-of- 
mass frame, this should be equal to the ionization potential 
of hydrogen, Jh = 13.6 eV. However, Soon adopts a thresh¬ 
old which is half of this value, for reasons which are unclear 


(see the lengthy discussion in Barklem 20071. In the high 
energy regime probed by experiments, this error is unim¬ 
portant, but at the low temperatures that we are interested 
in, the choice of threshold has a huge effect on the rate coeffi¬ 
cient, as we demonstrate below. In addition, |BarkIem] reports 
that the rate coefficient that one obtains by numerically in¬ 
tegrating the ionization cross-section given by Soon (using 
Soon’s choice of threshold) differs significantly at T < 10® K 
from the collisional ionization rate coefficient given by Soon. 

In view of these issues, we have chosen to proceed by 


ignoring the expressions given in Soon (19921 and instead 


analytically integrating the cross-section given in |Kunc fc| 
Soon (19911, using the correct value for the energy threshold. 
At the low energies of interest here, this yields a rate 


fcio,KS9i = 4.65 X IQ-^^T®''^ exp 


157800\ 3 -1 

—-— J cm s . (14) 


We compare these five different expressions for the H-H 
collisional ionization rate in Fig ure [4| We see that the rate 
derived from the Kune & Soon (19911 cross-section has an 


energy dependence that agrees well with the other versions 
of the rate coefficient. However, there is a significant dis¬ 
agreement when it comes to the normalization of the rate, 
with this value being around a factor of 100 smaller than the 
others at all temperatures. Fortunately, the impact of this 
uncertainty on Jerit is relatively small. Using a small value 
for the rate coefficient of reaction 10 has a similar effect to 
omitting it entirely, which as we saw in Paper I leads to a 
decrease in Jerit of around 20-30%. 

The remainder of the rate coefficients agree much bet¬ 
ter in terms of both their temperature dependence and their 
normalization. Even so, there remains a substantial scatter 
in the values at low temperatures. For example, as previ¬ 


ously noted, the Hollenbach & McKee (19891 and Lenzuni, 
Chernoff & Salpeter (19911 rates differ by around an order 


of magnitude at T = 7000 K, despite being based on the 
same atomic data. This scatter introduces a significant un¬ 
certainty into Jerit: in runs with a T4 spectrum, the value of 
Jerit varies by up to 80% depending on which version of the 
H-H collisional ionization rate coefficient we adopt, while 
in runs with a T5 spectrum, the uncertainty is even larger, 
around a factor of 2.5 



10000 

Temperature [K] 


Figure 4. Comparison of different rate coefficients for the col¬ 
lisional ionization of H by H (reaction 10). The values shown 
are taken from|Lenzuni, Chernoff &: Salpeter] (|1991[l (solid line), 
[Hollenbach fc McKee| ( |1989^ (dotte d line), [Hollenbach fc McKee| 
( [1979[ l (lower dashed iine),[D rawin ] l |1969[| (dash-dotte d line), and 
this work, based on the cross-section of [Kunc &: Soon[p991[ l (dot- 
dot-dot-dashed line). 


Finally, there remains the question of which of the dif¬ 
ferent rate coefficients one should choose for future calcula¬ 
tions. As the [Kune &: Soon] ( |1991| ) cross-section agrees very 
well with the available experimental data at E > 100 eV (see 
e.g. figure 3 in their paper), the rate coefficient based on this 
cross-section is probably the most accurate one available at 
the present time. However, the fact that there is a large gap 
between the range of energies for which the cross-section is 
experimentally constrained and the range of interest for de¬ 
riving the low-temperature form of the rate coefficient means 
that we must treat any of the possible rate coefficients with 
considerable caution. 


3.2.9 Collisional detachment of H by H (reaction 11) 


Considerable uncertainty remains in the rate of this reaction. 


An early study by Browne & Dalgarno (19691 gives two dif¬ 


ferent sets of values for the reaction rate, depending on what 
form is chosen for the interaction potential. If the [BardsIey,[ 
Herzenberg & Mandl ( 1966a|b I potential is adopted, then 
their numerical results are fit to within 10-20% by the ex¬ 
pression 


fcii.BD66 = 1.74 X 10-^"r°-^'’exp 


cm®s“\ (15) 


On the other hand, if the [Chen fc Peacher] ( [1968| ) potential 
is used, we instead obtain a rate coefficient that is well fit 
by the expression 

fell,CP68 = 1-74 X 10“^®T° ®'‘exp cm® s“\ (16) 

More recently, cross-sections for this reaction were given in 
the reviews of Janev, Langer & Evans (19871 and Janev, 


Reiter & Samm (20031. Based on the Janev, Langer & Evans 
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10 S. C. O. Glover 


(19871 data, Lenzuni, Chernoff & Salpeter (19911 propose a 
rate coefficient of the form 

fcii,LCS9i = 4.3 X cm^ s-\ (17) 


Abel et al. (19971 use the same cross-section data, but give 


a different analytical fit. At temperatures Te < 0.1, where 
Te = T/11604.4 K, their fit has the form 


7 o A \y 1A—Q/t-tI.78186 3 —1 

A:ii,A97 = 2.5do4 X 10 ie cm s , 


(18) 


while at higher temperatures, it has the form 


fcii,A97 = exp [—20.37260896 
+ 1.13944933(lnre) 

- 0.14210135(lnre)^ 

+ 8.4644554 x 10“^(lnre)® 

- 1.4327641 X 10“®(lnre)'‘ 

+ 2.0122503 X 10“'‘(lnre)® 

+ 8.6639632 x 10“®(lnre)® 

- 2.5850097 x 10“®(lnre)'^ 

+ 2.4555012 X 10“®(lnre)® 

-8.0683825 X 10"®(In Te)®] cm®s"\ (19) 

Finally, the more recent data from [Janev, Reiter fc Samm| 
( |2003[ ) can be used to derive a rate coefficient which is fit to 
within 1% over the temperature range 30 < T < 30000 K 
by the function 


fcii.jos = dex [—10.911243 

+ 2.0648286(logT3) 

- 0.40506856(logr3)^ 

+ 0.16462528(logT3)® 

- 0.10911689(logT3)'‘ 

+ 0.022537886(logr3)®l 


3 -1 
cm s , 


( 20 ) 


where T 3 = T/1000 K. 

In Figurej^ we show how these different versions of the 
rate coefficient behave as a function of temperature. We see 
immediately that at temperatures below a few thousand K, 
there is a large uncertainty in the value of the rate coeffi¬ 


cient. The values from Browne & Dalgarno (19691 fall off 


much more rapidly with decreasing temperature than the 
other versions of the rate coefficient. Moreover, even if we 
ignore this early work, and restrict our attention to the re¬ 
sults of more recent studies, we see that there is still a large 
uncertainty in the low temperature behaviour, with the rate 
coefficient that we have derived from the data in [Janev, | 
[Reiter fc Samm| ( |2003|) falling off more rapidly at low tem¬ 
peratures than the Lenzuni, Chernoff & Salpeter (19911 or 


Abel et al. (19971 rates. Fortunately, the chemical evolution 


of the gas is unlikely to be sensitive to this large uncertainty. 
The reason for this is that at the temperatures where there 
is an order of magnitude or more disagreement between the 
different rate coefficients, all of the expressions listed here 
predict values that are much smaller than the associative 
detachment rate. Therefore, in this regime, associative de¬ 
tachment dominates regardless of which expression we adopt 
for the rate of reaction 11 . 

Of greater concern is the smaller, but still substantial, 
disagreement between the values of the different expressions 
at temperatures close to 10'* K. In this regime, the differ- 
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Figure 5. Comparison of different rate coefficients for the colli- 
sional detachment of H“ by H (reaction 11). The values shown 
are taken from jAbel et al.j | |1997^ (s olid line), our fit tojjanev, Re-' 
[iter &: Sammj ( [2003^ (dotted li ne), [L enzuni, Chernoff &: Salpeter] 
( |1991| l (lower dashed line) and jBrowne &: Dalgarno| | |1969| ); in the 
latter case, the dash-dotted line shows the results that they obtain 
for the [Bardsley, Herzenberg &: Mandlj l |1966a|bt interaction po¬ 
tential, and the dash-dot-dot-dotted line shows their results for 
the [Chen &: Peach^ ( [1968^ potential. For comparison, we also 
plot the rate coefficient for reaction 3, the associative detachment 
of H“ by H (upper dashed line). The extremely large uncertainty 
in the rate of reaction 11 below a few thousand K is unlikely to 
significantly affect the chemical evolution of the gas, as in this 
temperature regime, all of the different expressions yield values 
that are much smaller than the associative detachment rate. 


ent determinations differ by a factor of three to fonr, and 
the magnitude of the rate coefficient is similar to that for 
associative detachment, meaning that this uncertainty po¬ 
tentially has a much larger impact on the H" abnndance 
and the H 2 formation rate. 

We have explored the effect that the uncertainty in this 
reaction has on the value of Jcrit. In runs with a T4 spec¬ 
trum, we find that the uncertainty introduced into our es¬ 
timate of Jcrit is very small, of the order of a few percent. 
This is unsurprising, as in these runs, the dominant destruc¬ 
tion process for H" is photodissociation and so even large 
changes in the rate of reaction 11 have little influence on the 
H~ abundance and hence little impact on the H 2 formation 
rate. In runs with a T5 spectrum, on the other hand, the 
uncertainty in the rate of reaction 11 has a much larger ef¬ 
fect, with Jcrit varying between 1630 (if we use the rate from 
Abel et al.|1997 1 and 2600 (if we use the rate from|Browne &: 


Dalgarno 1969 obtained with the Chen & Peacher potential). 


The total uncertainty in this case is therefore approximately 
60%. 


Although we would like to be able to recommend a 
particular choice of rate coefficient for this reaction, it is 
highly unclear which of the different possibilities is likely to 
be the most accurate. This process is therefore greatly in 
need of further theoretical or experimental study in order to 
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Simulating seed BH formation 11 


clarify the behaviour of its rate coefficient at temperatures 
T < 10‘‘ K. 


3.2.10 Collisional ionization of H by e (reaction 12) 

Several different expressions for the rate coefficient of this 
reaction are in reasonably common usage in the astrophysi- 
cal literature. Previous studies of the direct collapse model 
have typically adopted either the fit given in |Abel et al.| 
(19^ 


fcl2„ 


exp [-32.71396786 
+ 13.5365560(lnre) 

- 5.73932875(lnre)^ 

+ 1.56315498(lnre)^ 

- 0.28770560(lnre)'^ 

+ 3.48255977 x 10"^(lnTe)'^ 

- 2.63197617 x 10"^(lnTe)® 
+ 1.11954395 X 10"'‘(lnTe)'^ 
-2.03914985 x 10"®(lnTe)®l 


3 -1 
cm s , 


( 21 ) 


which is based on data from Janev, Langer & Evans (19871, 
or the expression given in Omukai ( 2001| , 

157800\ 


fci2,ooi = 4.25 X 10"“T^/^ exp - 




cm® s 


( 22 ) 


which comes from Lenzuni, Chernoff & Salpeter (19911. On 


the other hand, Hui & Gnedin (19971 give a different ex¬ 
pression 

21.11 ^_ A /2 


fcl2,l 


J'3/2 


(1 + [A/0.354]0'874)i.ioi 


cm® s ^(23) 


where A = 315614/T, which is their fit to Lotz (19671. Fi¬ 
nally, Scholz & Walters (19911 give the expression 


fcl2,SW91 = 


exp [-96.1443 
-f 37.9523 In r 

- 7.96885(lnr)® 

-f 8.83922 X 10"®(lnT) 

- 5.34513 X 10"®(lnT) 
-f 1.66344 X 10 

- 2.08888 X 10 


■®(lnT)® 
(InT)® 


-157800/r] cm®s"\ 


(24) 


which is based on the accurate cross-section measurements 


of Shah, Elliott & Gilbody (19871. 


The behaviour of these different expressions as a func¬ 
tion of temperature is compared in Figure At low tem¬ 
peratures (T <C 10® K), the collisional ionization rate has an 
exponential dependence on temperature, k \2 oc 
since only a small number of electrons in the exponential 
tail of the Maxwell-Boltzmann distribution have sufficient 
energy to ionize hydrogen. Therefore, in order to allow us 
to better see the difference between the different rate coeffi¬ 
cients, we do not plot the rate coefficient fei 2 in Figurej^ but 
instead plot the value of From the Figure, 

we see that at temperatures T > 10^ K, the different ver¬ 
sions of the collisional ionization rate coefficient agree fairly 
well. In this temperature regime, the different expression 
have roughly the same temperature dependence but disagree 
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Figure 6. Comparison of different rate coefficients for the colli¬ 
sional ionization of H by electrons (reaction 12). To highlight the 
difference between the various rate coefficients, we remove the 
underlying exponential temperature dependence by plotting the 
quantity where fci2 is the rate coefficient. The 


values shown in the plot come from 

Lenzuni, Chernoff Sz Salpeter 

(1991 

1 (solid line)|Abel et al.|<|1997 

\ (dotted line), 

Hui Sz Gnedin 

(1997 

1 (dashed line) and jScholz Sz Walters] ||1991 

1 (dash-dotted 


line). 


on the normalization of the rate coefficient by around 30%. 
At lower temperatures, the rates from [Lenzuni, Ghernoff fc 


[Salpeter] ( |1991[ ) , [Scholz &: Walters] ( |1991 1 and Hui & Gnedin 


(19971 continue to agree fairly well, while the Abel et al. 


(19971 rate begins to disagree strongly with the others. It 


seems likely that this is due to a problem with the |Abel[ 
[et al.[ fit at these low temperatures, although it should be 
remembered that this occurs in a regime where the value of 
fci 2 is falling off exponentially with decreasing temperature, 
and so for most applications, the [Abel et al.[ rate remains a 
suitable choice. 

We have determined the value of Jcrit that we obtain 
when using each of these four expressions for the rate coeffi¬ 
cient. In runs with a T4 spectrum, we find that the values all 
lie in the range Jcrit = 16.2-18.0, while in the case of a T5 
spectrum, the corresponding range of values is Jcrit = 1460- 
1630. The uncertainty in the rate of reaction 12 therefore 
introduces no more than a 10% uncertainty into our esti¬ 
mates of Jcrit. 

Finally, which of these rate coefficients should one ac¬ 
tually use in future calculations? Since the total uncertainty 


(statistical and systematic) in the Shah, Elliott & Gilbody 


(19871 measurements is unlikely to exceed 20% (see e.g. the 
discussion in jBarklem et al.|2011 1 and they are also in good 
agreement with available high-accuracy theoretical calcula¬ 
tions (e.g. Bray fc Stelbovics|1993 1, we recommend the use 
of the rate coefficient based on their measurements, i.e. the 
expression given in Scholz & Walters ( 1991|. 
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12 S. C. O. Glover 


3.2.11 Mutual neutralization of H with H'^ (reaction 22) 
In our fiducial model, we use the rate coefficient derived 


by Croft, Dickinson & Gadea (19991, based on the cross- 
sections of |Fussen Sz Kubacli| 1 1986p . This is well fit by the 
simple function 


fc22,C99 = 2.4 X ( 1 + 


20000) 


3 -1 
cm s 


(25) 


More recently, Urbain et al. (20121 state that their merged 
beam measurements “support the rate coeffcient compiled 
by Croft et ah”, but do not give numerical details of their 
measured cross-sections, meaning that it is not possible to 
assess how accurately their measurements agree with the 


earlier Fussen & Kubach (19861 cross-sections. The rate co¬ 


efficient for reaction 22 has also been calculated relatively 


recently by Stenrup, Larson & Blander (20091. Their results 


are well fit in the temperature range 10 < T < 10000 K by 
the function 

fc 22 .SLE 09 = 2.96 X 10"®T-^/^ - 1.73 x 10"® 


-h2.50 X 

- 7.77 X lO'^^Tcm® s“ 


(26) 


We compare this with the Croft et al. rate in Figure[7j where 
we plot the ratio A: 22 ,SLB 09 /fc 22 ,C 99 for temperatures in the 
range 10 < T < 10000 K. We see that the Stenrup et al. rate 
coefficient is higher than the Croft et al. rate coefficient at 
all temperatures, by around 20-25%. We have investigated 
the effect of using the Stenrup et al. rate coefficient in place 
of the Croft et al. rate coefficient in our models. We find that 
if we do so, then there is essentially no change in the value of 
Tcrit we recover for runs performed using a T4 spectrum. In 
runs with a T5 spectrum, however, we find that Jcrit = 1630 
when we use the Croft et al. rate coefficient and Jcrit = 1550 
when we use the Stenrup et al. rate coefficient. The current 
uncertainty in the value of the mutual neutralization rate 
coefficient therefore introduces an uncertainty of around 5% 
into our derived value of Jcrit. 

We note however that a number of treatments of pri¬ 
mordial chemistry in the literature use expressions for this 
rate coefficient based either on the measurements of|Moseley 


et al. {19701 or taken from the review of Dalgarno & Lepp 
(|1987p. These differ from the Croft et al. rate coefficient by 
factors of a few, as summarized in |Glover, Savin fc Jappsen| 
(20061. We have explored the effect on Jcrit of using one of 


these older values for the rate coefficient, and find that if we 
use the [Moseley et ah] value, we obtain Jcrit = 17.4 for a T4 
spectrum and Jcrit = 1070 for a T5 spectrum, while if we use 
the [Dalgarno fc Leppj value, we obtain instead Jcrit = 18.3 
for a T4 spectrum and Jcrit = 1940 for a T5 spectrum. We 
therefore see that the use of an older value for this rate coef¬ 
ficient can introduce an error of up to ~ 30% into the value 
we derive for Jcrit. 

Because of this, and in view of the fact that the most 
recent determinations of the rate coefficient for reaction 22 
agree well with the Croft et al. value and disagree signifi¬ 


cantly with the Moseley et al. {19701 and Dalgarno & Lepp 


(19871 values, we recommend that the rate coefficients for 


this reaction from the latter two sources should no longer 
be used when studying the direct collapse model. 



10 100 1000 lOOOC 

Temperature [K] 

Figure 7. Ratio of the rate coefficient for t he mutual neutral- 
ization of H“ by n't (reactio n 22) given by [Stenrup, Larson fcj 
[Elander1([2009[l to that given by[Croft, Dickinson fc Gadea[l[1999[|. 


4 IMPACT OF COOLING RATE 
UNCERTAINTIES ON Jcrit 

The cooling function that we use in our one-zone calcula¬ 
tions accounts for the effects of radiative cooling from a 
large number of different coolants, including atomic hydro¬ 
gen, atomic helium, H 2 , HJ, and HeH^ (see the detailed 
discussion in Glover &i Savin 20091. However, in practice. 


the only coolants that contribute substantially to the over¬ 
all cooling rate are atomic hydrogen and H 2 . We therefore 
focus our attention on these two coolants and do not con¬ 
sider the effects of uncertainties in the other cooling rates, as 
these introduce a negligible uncertainty into the total cool¬ 
ing rate. 


4.1 Ha cooling 

There are two main sources of uncertainty that we need to 
consider in the case of Ha cooling. At low densities, the cool¬ 
ing rate scales linearly with the collisional excitation rate, 
which is given by the sum of a series of contributions cor¬ 
responding to collisions between Ha and H, Ha and He, Ha 
and H’^, etc. In this regime, uncertainties in these contribu¬ 
tions directly affect the Ha cooling rate. At high densities, on 
the other hand, the Ha level populations reach LTE. In this 
regime, the Ha cooling rate is insensitive to the collisional 
excitation rate, and depends only on the partition function 
for the Ha molecule, the energies of the various excited ro¬ 
tational and vibrational levels, and their spontaneous radia¬ 
tive de-excitation rates. All of these quantities are known 
accurately, and so in principle there should be almost no 
uncertainty in the Ha cooling rate in these conditions. In 
practice, however, the Ha cooling rate in the LTE limit is 
often represented by a simple analytical function of tem¬ 
perature that approximates the true Ha cooling rate. Any 
error in this approximation therefore translates into an un- 
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Simulating seed BH formation 13 


Table 2. Sensitivity of Jcrit to the rate coefficients adopted for 
the listed cooling processes 


Process 

Sensitivity (T4) 

Sensitivity (T5) 

Lyman-0 

1.04 

1.54 

H2: collisions with H 

10.3 

4.53 

H2: collisions with He 

1.05 

1.03 

H2: collisions with H2 

1.00 

1.00 

H2: collisions with H"^ 

1.00 

1.00 

H2: collisions with e~ 

1.00 

1.00 


certainty in the resulting parameterization of the H 2 cooling 
rate. 


4-1.1 Low density limit 

In primordial gas, the most important collision partners re¬ 
sponsible for exciting the internal energy levels of the H 2 
molecule are H and He atoms, other H 2 molecules, protons 


and electrons (Glover & Abel 20081. In order to quantify 


which of these possible collisions partners is most important 
in the present case, we have carried out a similar analysis to 
that described in Section [sTl] We constructed a series of dif¬ 
ferent variants of our cooling model by adjusting the individ¬ 
ual collisional rate coefficients either upwards or downwards 
by a factor of lO'^ ® from their fiducial valuesj^We next used 
our variant models to determine Jcrit and then examined 
the difference in Jcrit that resulted from adjusting each rate 
coefficient either upwards or downwards. Finally, by taking 
the ratio of the largest and the smallest values of Jcrit that 
we obtained when modifying the rate of a particular colli¬ 
sion process, we could define a sensitivity for that process, 
analogous to the sensitivities determined for the different 
chemical rate coefficients in Section 

The results of our analysis are summarized in Table 
We see that Jcrit displays a high sensitivity to the rate of only 
one of our considered collisional excitation processes, the 
excitation of H 2 by collisions with H atoms. This strongly 
suggests that at the densities and temperatures where Jcrit 
is determined, collisions between H 2 and H dominate the 
total H 2 cooling rate. 

This behaviour differs from that found in metal-free pro¬ 
togalaxies that are not illuminated by extremely strong ra¬ 
diation fields, as in these systems H 2 -He collisions play an 
important role, as do H 2 -H^ and H 2 -e“ collisions if the gas 
is cooling from an initially ionized state (Glover & Abel 


2008). However, this difference is easy to understand. As we 


have already seen in paper I, in atomic cooling haloes illu¬ 
minated by a radiation field with J21 ~ Jcrit, the critical 
factor that determines whether or not the gas can cool is 
the behaviour of the H 2 fraction at densities n ~ 10 ® cm“® 
and temperatures T ~ 8000 K. In these conditions, the 
gas is mostly atomic and the fractional ionization is low, 
~ 10“^. If we compare the mean H 2 cooling rates per 
collision due to collisions with H, He, H"*" and e“ at this tem¬ 
perature and density, then we find that Anj.n/AHa.He ~ 4 


^ The latter were taken from |Glover &: Abel| | |200^ , although the 
rates for H2-H+ and H2-e“ collisions were updated to account 
for new data, as summarized in Appendix A of Paper I. 


and Ah 2 ,h/Ah 2 ,h+ ~ AHa.H/An^.e- ~ 0.3. Since the He:H 
ratio (by number) in primordial gas is approximately 0.08, it 
is clear that collisions with helium contribute only a few per¬ 
cent of the total H 2 cooling rate, and that the contributions 
from collisions with protons and electrons are negligible. In 


the scenario considered by Glover & Abel (20081, gas cools 


and recombines from an initially ionized state, but there is 
no ultraviolet background. In this case, H 2 cooling becomes 
effective much earlier in the collapse, when the fractional 
ionization is significantly larger. In addition, [Glover fc Abel| 
(20081 were primarily concerned with the behaviour of the 


cooling rate at low temperatures, where the difference be¬ 
tween the contribution from H 2 -H collisions and from colli¬ 
sions with He, H'^ or electrons is much more pronounced. 

Since our analysis demonstrates that it is essentially 
just the H 2 -H collisional excitation rate that plays a role in 
determining Jcrit, we choose to focus our attention on this 
process, and will not discuss any uncertainties that may exist 
in the other collisional excitation rates. In our default model, 
the treatment we use for the H 2 -H contribution to the H 2 


cooling rate is taken from Glover & Abel (2008). They give 


the following analytic fit to the H 2 cooling rate, valid in the 
temperature range 1000 < T < 6000 K: 


A; 


H2,GA08 = 


dex [-24.311209 
+ 4.6450521(logT3) 

- 3.7209846(logT3)® 
-b 5.9369081(logT3)® 

- 5.5108047(logT3)'‘ 
-b 1.5538288(logT3)®l 


erg s 


(27) 


where T 3 = T/IOOO k[^ This fit is based on the set of colli¬ 
sional excitation rate coefficients computed byjWrathmall fcj 
Flower (2007; see also |Wrathmall, Gusdorf fc Flower|2007 1. 
These authors use a fully quantal treatment of the H 2 -H 
system in their calculations and make use of the potential 


energy surface of Mielke, Garrett & Peterson (2002). 


The other main treatment of H 2 cooling that is in 
widespread use in numerical models of primordial gas is that 


of Galli & Palla (19981. They give the following analytical 


fit to the H 2 -H cooling rate: 

Ah 2 ,gp 98 = dex [-103.0-b 97.59logr- 48.05(logr)^ 
- 48.05(logr)^ -b 10.80(logr)® 


— 0.9032(logr)"^l ergs ® cm®. 


(28) 


This fit is based on the quantal calculations of jForrey et al.j 
(19971 at T < 600 K and the quasi-classical calculations of 


Mandy & Martin (19931 at T > 600 K. In both cases, older 


versions of the H 2 -H potential energy surface were used, as 
summarized in Galli fc Palla[ (1998 1 . 

Finally, an alternative version of the H 2 -H cooling rate 
is given in Martin, Schwarz & Mandy (19961. These authors 


present a complicated fit to the H 2 cooling rate in the regime 
where H collisions dominate that is valid over a wide range 
of different temperatures. However, the low density limit of 


Glover & Abel 


12008 I also give fits for in the temperature 


ranges 10 < T < 100 K and 100 < T < 1000 K, but the behaviour 
of the H2 cooling function at these temperatures has no influence 
on the value of Jcrit. 
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14 S. C. O. Glover 


this fit is easy to extract and can be written as 


Ah 2 ,msm 96 = dex [+195.09 — 137.9986 logT 
+ 29.802(logr)^ 

- 2.165163(logr)® 

+ 91.60106 log(1.0 + 2680.075/r) 

- 105.7438 log(1.0 + 2889.762/r) 

-4328.147/r] ergs“^cm^ (29) 


This version of the cooling rate is also based on the calcula¬ 
tions of Mandy & Martin ( 1993[ ). 

In the upper panel of Figure we compare the tem¬ 
perature dependence of these three versions of the H 2 -H 
cooling rate. We see that although the different rates dis¬ 
agree somewhat at temperatures close to 1000 K, there is 
good agreement in the temperature range T ~ 6000-8000 K 
relevant for Jcrit- At even higher temperatures, T > 10'* K, 
there is significant disagreement between the rates, but this 
is simply a consequence of the breakdown of the [Glover ^ 
|Abel| (|2008|) and |Galli & Palla| (|1998|) analytical fits in this 
temperature regime; as previously noted, the [Glover fc Ab'el 


fit is formally only valid at T < 6000 K, while the Galli & 
Palla fit is valid only for T < 10* K. Fortunately, at these 


very high temperatures, the precise behaviour of Aua.u is 
unlikely to be important, as Lyman-a cooling dominates in 
this temperature regime. 

As the H 2 cooling rate varies over several orders of mag¬ 
nitude in the temperature range that we are examining, it is 
difficult to tell from the upper panel of Figure [^exactly how 
well the different cooling rates agree at T ~ 6000-8000 K. 
Therefore, in the lower panel of the Figure, we plot instead 
the ratio of the lGalli fc Palla! rate to the lGlover fc Ab^ rate 
(solid line) and the ratio of the [Martin, Schwarz fc Mandy 


rate to the Glover & Abel rate (dashed line). We see that at 
the temperatures of interest, all three rates agree to within 
around 10 %. 

We have investigated the impact of this small uncer¬ 
tainty on Jcrit. We find that in runs with a T4 spectrum, 
it introduces an uncertainty of roughly 10 % into our deter¬ 
mination of Jcrit, while in runs with a T5 spectrum, the 
corresponding uncertainty is ~ 5%. 


4-1.2 LTE limit 


In the LTE limit, the H 2 cooling rate can be computed very 
accurately, since it only depends on the energies and ra¬ 
diative de-excitation rates of the rotational and vibrational 
levels of H 2 , all of which are known with a high degree of pre¬ 
cision. Glover (2011 1 numerically evaluated the LTE cooling 


rate and produced the following analytical fit for the cool¬ 
ing rate coefficient, which is valid in the temperature range 
100 < T < 10000 K and has an error of no more than 5% 
over this range: 


Ah 2 ,ltb,gii = dex [—20.584225 

+ 5.0194035 log T 3 

- 1.5738805(logr3)^ 

- 4.7155769(logr3)^ 
+ 2.4714161(logr3)* 
+ 5.4710750(logr3)® 




1000 10000 
Temperature [K] 


Figure 8. Upper panel: comparison of H2-H cooling rate coeffi¬ 
cients.The solid line shows the high temperature fit from [Glover 
Sz, Abelj l |2008[| , the dotted line shows the expression from |GaIIi 
Sz Palla] (|1998|l and the da shed line shows the expression from 
Martin, Schwarz &: Mandy] p'99^ . The significant disagreement 
in behaviour at T > 10^ K is a consequence of the fact that the 
fits given in [Glover &: Abel| l |2008| | and [Galli Sz Pallaj < [1998|l are 
not valid at these temperatures. Lower panel: ratio of the [Galli[ 
[&: Palla] | |1998[| rate to the[Glover &: Abel[([2008| l rate (solid line) 
and of t he [Martin, Schwarz Sz Mandy[ < [1996| l rate to the [Glover[ 
\Sz Abel[ ||2008[| rate (dashed line). 


-3.9467356(logr3)® 

- 2.2148338(logr3)'^ 

+ 1.8161874(logr3)®] ergs“\ (30) 


However, many treatments of H 2 cooling in primordial gas 
(see e.g. Galli fc Palla||19M Bryan et al.|2014 Grassi et al. 


20141 continue to make use of the approximate analytical 


expression for the H 2 LTE cooling rate given in [Hollenbach[ 


& McKee (19791: 


A 


H2,LTB,HM79 = 


— Ah 2 .rot + Aho 


(31) 
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where 

+ 3 X exp(-0.51/r3) erg s"\ (32) 

and 

AH 2 ,vib = 6.7 X 10 '^®exp(- 5 . 86 /T 3 ) 

+ 1.6 X 10"^®exp(-11.7/r3) ergs"\ (33) 


We compare these two cooling rates in Figure where 
we plot their ratio, Ah 2 ,ltb,gii/Ah 2 .lte,hm 79 . We see that 
at temperatures below T ~ 2000 K, the two rates agree 
reasonably well, with differences of no more than 10%. At 
higher temperatures, however, the expression from |Hollen-| 
bach & McKee (19791 yields a larger cooling rate than the 
one from [Glover ( |201ip , with the discrepancy worsening to 
as much as a factor of two as we near 10^ K. 

In practice, this discrepancy does not have a large im¬ 
pact on our estimate of Jcrit, as this is determined at densi¬ 
ties somewhat lower than the H 2 critical density. If we com¬ 
pare the results of runs performed using the [Glover (|2011 l 
LTE rate with those performed using the|Hollenbach fc Mc-| 


Kee (19791 rate, we find that Jcrit differs by around 10%. 


This holds regardless of whether we use a T4 or a T5 spec¬ 
trum. 


4.2 Lyman-a cooling 

4-2.1 Collisions with electrons 


Most numerical models of primordial gas use the following 
expression for the cooling rate coefficient arising from H-e“ 
collisions (commonly referred to simply as Lyman-a cool¬ 
ing): 


A 


H,C92 


7.5 X 10"^® 


1 + T, 


1/2 


-118348/T -1 3 

e ' erg s cm 


(34) 



Figure 10. Comparison of the Lyman-a cooling rate coeffi¬ 
cients from [Cen| l |1992] | (solid line) and |Scholz &: Walters] l |l991| | 
(dashed line). Both rate coefficients have been divided by a fac¬ 
tor in order to remove the underlying exponential 

dependence on temperature, allowing us to more clearly see the 
difference in behaviour between the two expressions. 


where T 5 = T/10^ K. This expression comes from Cen 


(19921, and is a minor modification of the cooling rate coef¬ 


ficient originally given in Black (19811. In order to assess the 
level of uncertainty in this parameterization of the cooling 
rate, we compare it with the expression given in |Scholz «S| 
Walters! (|1991|): 


Ah,sw 91 = 10 ^“/(T)e ergs ^ cm^, 

where 


(35) 


/(T) = exp [213.7913- 113.9492 InT 

+ 25.06062(lnr)^ 

- 2.762755(lnr)^ 

+ 0.1515352(lnr)'‘ 

- 3.290382 X lO'^lnT)®] . (36) 

These two versions of the Lyman-a cooling rate are com¬ 
pared in Figure [To] To enable us to more clearly distinguish 
the differences in the behaviour of the two rates, we have 
divided both of them by a factor g-ii*348/T order to take 
out the underlying exponential dependence on temperature. 
We see that although there is some disagreement in the rates 
at low temperatures (where Lyman-a cooling is unimpor¬ 
tant), at the temperatures relevant for Jcrit these two differ¬ 
ent prescriptions agree to within a few percent. This small 
difference introduces an uncertainty of around 1 % into our 
estimate of Jcrit. 


4-2.2 Collisions with H and He atoms 

Electronic states of atomic hydrogen with n > 2 can also be 
excited by collisions with H and He atoms. The rate coef¬ 
ficients for these processes are much smaller than for H-e“ 
collisions, but as in the case of the collisional ionization of 
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hydrogen, it is conceivable that these processes may become 
important if the fractional ionization of the gas is very small. 

At the temperatures of interest in this stndy, the atomic 
hydrogen cooling rate is dominated by the excitation of the 
2 p and 2s levels, and so there is no need for ns to consider 
excitation to states with n 3. We can therefore write the 
cooling rate coefficient that represents the effect of H-H col¬ 
lisions as 


Ah,h = 10.2eV X gi 2 , ergs ^cm^, 


(37) 


where qi 2 is the collisional excitation rate coefficient for tran¬ 
sitions from n = 1 to n = 20 

As in the case of H-H collisional ionization, there are 
very few studies of H-H collisional excitation that consider 
the low energies relevant for our present study. One possibil¬ 
ity is the rate given in |Lenzuni, Chernoff fc Salpeter (19911, 
which is based on the work of Drawin ( 1969[ ): 

qi2 = 5.8 X Q ^ gg ^ 10"®T) 


, 118348\ 3 -1 

X exp (- - — 1 cm s , 




(38) 


and so 


Ah,h,lcs 9 i = 9.5 X 10-2®T^''^ (l.O -f 1.69 x 10"®t) 

/ 118348\ _i 3 
X exp I-—— I erg s cm . (39) 

Collisional excitation of H by H is also treated by|Soon| 


119921. However, in this case we encounter the same problem 


as we did when examining the treatment of H-H collisional 
ionization by the same author: the energy threshold adopted 
in Soon’s calculations is only half of the value that it should 


be (Barklem 20071. If we correct for this problem and nu¬ 


merically integrate the cross-section given in Soon (19921 


using the correct energy threshold, then the cooling rate co¬ 
efficient that we obtain can be fit to within 10% over the 
temperature range 1000 < T < 20000 K by the function 

118348 \ 


Ah,H,S 92 b = 5.8x10 °'^®exp 


T J 


erg s 


(40) 


In Figure pH| we compare these two different H-H cool¬ 
ing rate coefficients. At the temperatures relevant for Jcrit, 


the rate coefficient from Lenzuni, Chernoff & Salpeter (19911 


predicts about an order of magnitude more cooling than the 
rate coefficient based on Soon| (19921. If we compare these 
values to the H-e“ cooling rate coefficients discu ssed ear¬ 
lier, then we find that Ah h/Ah^b- ~ 10 ® for the 
Chernoff & Salpeter (19911 rate coefficient and ~ 


Lenzuni, 
® for 


10 “ 


the rate coefficient based on Soon (1992 l. As the actual frac¬ 


tional ionization of the gas in the conditions where Jcrit is 
determined is around a; ~ 4 x 10“® (see e.g. Paper I, Fig¬ 
ure 1), it is unsurprising that if we include the effect of H-H 
cooling using either of the two rate coefficients from the lit¬ 
erature, we find that the effect on Jcrit is negligible. The 
large uncertainty in the H-H cooling rate therefore has no 
influence on the value of Jcrit- 

As far as the excitation of H by collisions with He atoms 
is concerned, Lenzuni, Chernoff & Salpeter ( 1991| give the 


® In principle, one should distinguish between the Is —> 2s and 
Is —t 2p transitions, but in practice, treatments of this process in 
the literature often assume that both transitions share the same 
rate coefficient. 
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Figure 11. Comparison of the H-H collisional excitation cooling 
rate coefficients from [Lenzuni, Chernoff fc Salpeter] ( |1991| (solid 
line) and [Soon] ||1992| (dashed line). 


following expressions for the collisional excitation rates of 
the Is —t 2 p and Is —>■ 2s transitions respectively: 

gi^2p = 1.2 X exp cm®s“\(41) 

ql-^2s = 1.2 X 10~^^r^'^ exp cm® s“\ (42) 


These collisional excitation rates are based on the authors’ 
fits to the data of Grosser & Kruger (19841. Using these ex¬ 
citation rates, we can easily derive the corresponding cooling 
rate coefficient: 

196250\ 


Ah,H e = 2.2 X 10“®®r®-^exp (- 




ergs ® cm®. (43) 


At T ~ 10'* K, this is approximately 10® times smaller than 
the expressions for Ah given above. Therefore, the cooling 
rate due to H-He collisions becomes comparable to that due 
to H-e“ collisions only when the fractional ionization of the 
gas is very small, x ~ 10“®. In the conditions relevant for 
Jcrit, the fractional ionization is around a factor of 40 higher 
than this, and so it is safe to conclude that H-He collisions 
make a negligible contribution to the total cooling rate. 


5 DISCUSSION 

5.1 Combining the uncertainties 


It is useful at this point to draw a distinction between the 
uncertainties in Jcrit that exist due to the simplifications 
made in our chemical model and those that exist due to un¬ 
certainties in the chemical rate coefficient data. The case 
of H 2 self-shielding provides a good example of the for¬ 
mer type of uncertainty. As we have already discussed, the 
value of Jcrit that we recover in runs with a T5 spectrum is 
highly sensitive to our choice of H 2 self-shielding function, 
increasing by more than a factor of five if we use the origi¬ 
nal Draine & Bertoldi (19961 self-shielding function instead 
of the more recent version given in [Wolcott-Green, Haimml 
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Table 3. Error in Jcrit introduced by uncertainty in the rate 
coefficients of the listed reactions 


No. 

Reaction 

Uncertainty 
T4 spec. 

in Jcrit (%) 
T5 spec. 

2 

H2-|-H->H-|-H-|-H 

35 

25 

3 

H- -h H ^ H2 -f e- 

40 

40 

6 

H + e“ ^ H“ + 7 

20 

30 

10 

H-|-H-5-H+-|-e--|-H 

80 

120 

11 

H--|-H-s>H-|-H + e- 

Negligible 

60 


& Bryan (20111. However, this large uncertainty ultimately 


stems from the fact that these self-shielding functions are 
both approximations, valid in different limits. In reality, it 
is likely that neither alternative provides a completely accu¬ 
rate description of the behaviour of the H 2 in the gravita¬ 
tionally collapsing gas over the whole range of number densi¬ 
ties of interest. Since the molecular data needed to compute 
the H 2 photodissociation rate and self-shielding function is 
known accurately, we can in principle eliminate this uncer¬ 
tainty simply by using a more sophisticated treatment of 
the H 2 in our simulation that tracks the population of the 
different rotational and vibrational levels and computes the 
photodissociation rate separately for each level. Similarly, 
any uncertainties in the H 2 photodissociation rate that are 
introduced by the simple approximations commonly used to 
compute N 112 can be eliminated by abandoning these ap¬ 
proximations using a more accurate approach to compute 


the H 2 column density (see e.g. |Wolcott-Green, Haiman fc 


Bryan|2011 Hartwig et al.|[^015i. 

Uncertainties in Jcrit that exist due to uncertainties in 
the chemical rate coefficient data are less easily dealt with. 
In order to reduce their impact, it is necessary to determine 
the rate coefficients more accurately, either experimentally 
or computationally, but this can be an extremely challeng¬ 
ing and time-consuming occupation. It is therefore useful to 
assess the size of the uncertainty in Jcrit that results from 
these rate coefficient uncertainties, so that we can better 
assess how important it is to work on reducing these. 

From the discussions in Sections]^ andit is clear that 
in runs performed with a T4 spectrum, there are only four 
chemical reactions that introduce significant uncertainties 
into Jcrit: the collisional dissociation of H 2 (reaction 2), the 
associative detachment of to form H 2 (reaction 3), the 
formation of H“ by radiative association (reaction 6) and 
the collisional ionization of H by H (reaction 10). These re¬ 
actions individually introduce uncertainties of at least 25% 
into Jcrit, as summarized in Tablewhile the uncertainties 
contributed by the other processes discussed in Sections 
and 1^ are at the level of 10% or less. In runs with a T5 
spectrum, the uncertainties introduced by these four reac¬ 
tions remain important, but Jcrit also depends strongly on 
the rate of reaction 11, the collisional detachment of H“ by 
collisions with H. 

We can place some bounds on the total uncertainty in¬ 
troduced into Jcrit by these reactions by deliberately choos¬ 
ing rate coefficients from amongst the different possibilities 
that either maximize or minimize Jcrit. If we do this, then 
we find that in runs with a T4 spectrum, Jcrit,min = 9.8 and 
Jcrit,max = 47.5. In runs with a T5 spectrum, we find instead 
that Jcrit,min = 850 and Jcrit,max = 5540. Therefore, in the 


worst case, chemical rate coefficient uncertainties could in¬ 
troduce an uncertainty of around a factor of five into our 
estimates of Jcrit. 

In practice, the situation is probably not as bad as this, 
as the uncertainties in the different rate coefficients will can¬ 
cel with each other to some extent. To assess the impor¬ 
tance of this, we recalculated Jcrit a dozen times for both 
the T4 and T5 cases, in each calculation using rate coeffi¬ 
cients for reactions 2, 3, 6, 10 and 11 that were randomly 
selected from amongst the different possibilitiesl^This pro¬ 
cedure gave us 12 different values of Jcrit for each choice of 
spectrum, and we then calculated the mean and standard 
deviation of these sets of values. We found that for runs 
with a T4 spectrum, Jcrit = 21 ± 5, while for runs with a 
T5 spectrum, Jcrit = 2400 ± 700. We can therefore conclude 
that the total uncertainty in Jcrit due to uncertainties in the 
chemical rate coefficients probably does not exceed a factor 
of two and certainly does not exceed a factor of five. 

For comparison, determinations of Jcrit made using 3D 
simulations rather than the simplified one-zone approach 
used here typically find that Jcrit varies by at least a fac¬ 
tor of a few from halo to halo owing to minor differences 
in the details of the collapse of the gas and its resulting 
temperature structure (Shang, Bryan & Haiman|2010 Latif 
eTaLllMI l. Even larger uncertainties are introduced if we 
abandon the use of our simplified T4 and T5 spectra and 
consider more realistic models for the interstellar radiation 


field I Sugimura, Omukai fc Inoue|2014 Agarwal & Khochfar 


2015 Agarwal et al.|2015 l, particularly once we account for 

the fact that there may be a non-negligible X-ray compo¬ 
nent ( Inayoshi fc Tanaka|2015' Latif et al.|2015 (. Therefore, 
it seems clear from the results of our study that astrophysi- 
cal uncertainties have a larger impact on Jcrit than chemical 
uncertainties. 


5.2 Comparison with previous work 


Previous studies of the direct collapse model using a one- 
zone approach have recovered a range of different values for 
Jcrit, with values ranging from 18-39 in runs with a T4 spec¬ 
trum, and 1400-16000 in runs with a T5 spectrum (see the 
overviews of previous work in [Inayoshi fc Tanaka|[^15| and 
Agarwal et al.|[2015 |. Some of this uncertainty, particularly 
in the T5 runs, is due to differences in the treatment of 
H 2 self-shielding ( [Sugimura, Omukai fc Inoue|[^14p , and 
some can be ascribed to differences in the composition of 
the chemical networks used to model the gas, as explored 
in Paper I. However, it is also useful to assess how much 
of this uncertainty might reasonably be due to differences 
in the chemical rate coefficients used within these different 
models. 

Specifically, in this section we compare the results that 
we obtain using our standard choices for the various rate co¬ 
efficients with those that we obtain using the set of rate co¬ 
efficients adopted within (a) the enzo primordial chemistry 


In the case of reactions 2 and 3, we accounted for the 25% 
systematic uncertainty in the two rate coefficients by considering 
both our fiducial versions and also variants that were 1.25 and 
0.75 times these fiducial versions. 
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Table 4. Values of Jcrit derived using the ENZO and KROME pri- Jcrit for the T4 spectrum comes from differences in 

mordial chemistry networks the chemical models used for the different studies. 




•^crit 

Model 

T4 spec 

T5 spec. 

ENZO (unmodified) 

34.0 

3050 

ENZO (modified) 

24.1 

2480 

KROME (unmodified) 

8.2 

690 

KROME (modified) 

11.2 

1080 


network ( Bryan et al.|2014|), and (b) the krome astrochem- 
istry package ( Grassi et al.|2014 1 . We choose to compare our 
results with these two alternative treatments because many 
of the previous studies of direct collapse use one or the other 
of them, and moreover it is easy to gather information on 
exactly which rate coefficients are used within these models. 


5.2.1 The ENZO network 

Direct comparison between the results we obtain with our 
chemical network and those obtained using the ENZO primor¬ 
dial chemistry network is complicated by the fact that the 
ENZO network not only makes different choices for a num¬ 
ber of the chemical and cooling rate coefficients, but also 
consists of a different set of chemical reactions from those 
included in our model. In particular, the chemical network 
implemented in the most recent version of the ENZO code 
(version 2.4) omits two processes that have a significant in¬ 
fluence on the value of Jcrit: the collisional ionization of H by 
H (reaction 10) and the dissociative tunneling contribution 
to the H 2 collisional dissociation rate. Therefore, to allow 
us to distinguish between the differences in Jcrit caused by 
the omission of these reactions and the differences in Jcrit 
caused by the differences in the rate coefficients, we compare 
our results with those from two additional sets of calcula¬ 
tions. In one, we use the unmodified ENZO network. In the 
other, we use a version of this network that we have modified 
to include the two omitted processes mentioned above. 

The results of our comparison are shown in Table We 
see that the enzo network systematically produces higher 
values of Jcrit that those derived using our network. With 
the unmodified ENZO network, Jcrit is increased by almost a 
factor of two for both the T4 and the T5 spectrum. With 
the modified network, on the other hand, the increase is 
smaller: around 30% for the T4 spectrum and 50% for the 
T5 spectrum. Note that these values were computing using 
the case A hydrogen recombination rate, which is the default 
choice in ENZO. If we instead use the case B rate, which is 
more physically appropriate in the present case, we recover 
values of Jcrit that are almost a factor of two larger than 

the literature for Jcrit in the 
zone calculation comes from 
), who find that Jcrit = 39, 
roughly a factor of two larger than most other estimates. 
However, they use the enzo chemistry network in their 
study, and our results here suggest that much of this factor 
of two differences is due to the composition of their network 
and their choice of reaction rate coefficients. It therefore 
seems plausible that most of the scatter in the reported val- 


those in Tabled 

The highest value quoted in 
case of a T4 spectrum and a one- 


Shang, Bryan & Haiman (2010 


5.2.2 The krome package 

As already noted in Paper I, the krome package provides 
a number of different networks suitable for simulating pri¬ 
mordial gas, but provides little guidance regarding which 
network should be used for which application. For studying 
direct collapse, the most suitable choice appears to be the 
react^xrays network. We therefore use this as the basis of 
our comparison]^ As in the case of the enzo network, this 
omits the effects of the collisional ionization of H by H, and 
so when carrying out our comparison, we consider both the 
unmodified version of the react-xrays network plus a modi¬ 
fied version that includes this reaction. 

Our results are shown in Table ID We see that the un¬ 
modified network underestimates Jcrit by around a factor of 
two. Including the collisional ionization of H by H improves 
the situation slightly, but even in this case the krome net¬ 
work still significantly underestimates Jcrit- We have inves¬ 
tigated the causes of this discrepancy and find that it is 
largely driven by the choice of recombination rate. The 
krome react-xrays network adopts the case A rate for this 
process, which causes the H"*" abundance to fall off more 
rapidly with increasing density than in our models. If we 
instead use the more physically appropriate case B rate, we 
find much better agreement between the krome results and 
our own determination of Jcrit- This demonstrates that the 
error introduced into Jcrit is solely a consequence of the par¬ 
ticular choice of network and rate coefficients, and not any 
intrinsic problem within the krome infrastructure itself - if 
KROME is used simply as an ODE solver and cooling func¬ 
tion, with a user-supplied chemical network that includes all 
of the important chemical reactions and an up-to-date set 
of rate coefficients, then we expect it to be able to produce 
exactly the same values of Jcrit as we find in our study. 


6 SUMMARY 


In this paper, we have attempted to determine the extent to 
which estimates of Jcrit are affected by uncertainties in the 
chemical rate coefficients and cooling rates used to model 
the chemical and thermal evolution of primordial gas. To 
do this, we first examined how sensitive Jcrit is to varia¬ 
tion in the rate coefficients of the 26 reactions included in 


the reduced chemical network of Glover (20151. For those 


reactions where the sensitivity is large, we then critically 
examined the information available regarding the values of 
their rate coefficients, and quantified the effect of uncertain¬ 
ties in these values on our derived value of Jcut- We also 
carried out a similar analysis for the main cooling processes 
included in our model. 

We found that there are five key chemical reactions 
where the combination of sensitivity and uncertainty leads 
to an appreciable (> 20%) uncertainty in Jcut- These are 


® Specifically, we consider the version present in the main KROME 
repository on March 26th, 2015, at which time the most recent 
commit was cf 3d6b882d35545657d0a5daa7965b5c62c7e970 
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the collisional dissociation of H 2 by H (reaction 2), the for¬ 
mation of H 2 by associative detachment of H~ (reaction 3), 
the formation of H~ by radiative association (reaction 6), 
the collisional ionization of H by H (reaction 10) and the 
collisional detachment of H“ by H (reaction 11), although 
the last of these is important only in situations where H“ 
photodetachment is unimportant. Amongst these reactions, 
the single largest uncertainty comes from the collisional ion¬ 
ization of H by H, a process which is very poorly constrained 
at the low energies relevant for this study | Barklem 20071. In 
addition, in gas illuminated with a hard T5 spectrum, Jcrit 
is also highly sensitive to the value of the H 2 photodissoci¬ 
ation rate, although in this case the intrinsic uncertainty in 
the rate is small, and the main difficulty comes from the fact 
that computationally efficient treatments of H 2 self-shielding 


are often highly over-simplihed (see e.g. Wolcott-Green & 


Haiman 20111. The total uncertainty introduced into Jcut 


by uncertainties in the remaining 20 reactions is small in 
comparison to the effects of these six reactions. Similarly, 
the remaining uncertainties in the cooling rates of H 2 and 
atomic hydrogen do not signihcantly affect Jcrit- 

In the unlikely event that the errors in the different 
rate coefficients combine so as to maximum or minimize 
Tcrit, we find that for a T4 spectrum, Jcrit,min = 9.8 and 
Jcrit,max = 47.5, while for a T5 spectrum, Jcrit,min = 850 
and Jcrit,max = 5540. Therefore, in the worst case, chemical 
rate coefficient uncertainties could introduce an uncertainty 
of around a factor of hve into our estimates of Jcrit. In the 
more likely case that the errors in the different rate coeffi¬ 
cients are uncorrelated, we find that Jcrit probably lies in the 
range Jcrit = 21±5 for a T4 spectrum and Jcrit = 2400±700 
for a T5 spectrum. 

Finally, we have compared the values of Jcrit that we 
have derived using our reaction network and rate coefficients 
with those that we obtain if we use instead the network 
and rate coefficients from (a) the enzo hydrodynamical code 
or (b) the krome astrochemistry package, both of which 
have been used in a number of previous studies of the direct 
collapse model. We find that the enzo chemical model tends 
to overestimate Jcrit by around a factor of two, although this 
discrepancy increases to a factor of four if the physically 
appropriate case B hydrogen recombination rate is used in 
place of the case A rate that is the default in enzo. The 
krome model, on the other hand, underestimates Jcrit by 
around a factor of two. 
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